Determining bandwidth parameter for RKHS

For RKHS, h is bandwidth parameter. NOT heritibility.

Based on Campos et al., 2010,

“As shown in Fig. 4, the value of the bandwidth parameter that gave the best predictive ability was in the range (2,4), except for environment E2 in which values of h near one performed slightly better.”

where figure 4 is the Estimated MSE (mean-squared error) between CV predictions (yHat) and observations (y) versus values of the bandwidth parameter, h, from 0.1 - 10. The bandwidth parameter with the lowest MSE gave the best predictive ability.
So as we look the the PMSE plot against h, we should select the bandwidth parameter h that gives the lowest MSE for each of your pheno datasets.

I need to determine what is the optimal bandwidth parameter to use for each dataset. Campos et al., 2010 should shed some light on how to select the h.

## Determining bandwidth parameter for RKHS ----------------
setwd("D:/PHS Genomic Selection Project/Data Analysis/GS/20181218")
# setwd("/Users/shantel/Desktop/20181218 GS Mac")
folds=5              
run=5
nIter = 12000
burnIn=2000

#Bandwidth parameters to compare
h=c(.01,.1,.4,.8,1.5,3,5)

#Your Data ~~~~~~~~~~~~~~~~
myY = myYc; myGD = myGDc
y = myY[,2] 
X = myGD; n<-nrow(X); p<-ncol(X)
D<-(as.matrix(dist(X,method='euclidean'))^2)/p

#Statistic Parameters
PMSE<-numeric() ; VARE<-numeric(); VARU<-numeric() ;
pD<-numeric(); DIC<-numeric()
fmList<-list()

#RKHS Calculation ~~~~~~~~
for(i in 1:length(h)){
  K<-exp(-h[i]*D)
  ETA<-list(list(K=K,model='RKHS'))
  prefix<- paste(h[i], "_",sep="")
  fm<-BGLR(y=yNA,ETA=ETA, nIter=nIter,burnIn=burnIn,saveAt=prefix)
  fmList[[i]]<-fm
  PMSE[i]<-mean((y[tst]-fm$yHat[tst])^2)
  VARE[i]<-fm$varE
  VARU[i]<-fm$ETA[[1]]$varU
  DIC[i]<-fm$fit$DIC
  pD[i]<-fm$fit$pD
}

#Plots for each h ~~~~~~~
R2<-1-PMSE/mean((y[tst]-mean(y[-tst]))^2)

plot(VARE~h,xlab="Bandwidth", 
     ylab="Residual Variance",type="o",col=4)
plot(I(VARE/VARU)~h,xlab="Bandwidth",
     ylab="variance ratio (noise/signal)",type="o",col=4)
plot(pD~h,xlab="Bandwidth", ylab="pD",type="o",col=2)
plot(DIC~h,xlab="Bandwidth", ylab="DIC",type="o",col=2)
plot(R2~h,xlab="Bandwidth", ylab="R-squared",type="o",col=2)
plot(PMSE~h,xlab="Bandwidth", ylab="PMSE",type="o",col=2)

Also, I just found another article discussing selection of the bandwidth parameter Pérez-Elizalde et al., 2015. Need to review, but for now, I think we have a good start for h.

LS0tDQp0aXRsZTogImg6IEJhbmR3aWR0aCBQYXJhbWV0ZXIiDQpvdXRwdXQ6IA0KICBodG1sX2RvY3VtZW50Og0KICAgIGNvZGVfZG93bmxvYWQ6IHllcw0KICAgIGRmX3ByaW50OiBwYWdlZA0KICAgIHRoZW1lOiBmbGF0bHkNCi0tLQ0KDQpbRGV0ZXJtaW5pbmcgYmFuZHdpZHRoIHBhcmFtZXRlciBmb3IgUktIU10oaHR0cHM6Ly9naXRodWIuY29tL3NoYW50ZWwtbWFydGluZXovQ29ybmVsbF9QSFMvYmxvYi9tYXN0ZXIvUEhTJTIwR2Vub21pYyUyMFNlbGVjdGlvbiUyMFByb2plY3QvRGF0YSUyMEFuYWx5c2lzL0dTLzIwMTgxMjE3L1BIUyUyMEdTJTIwdjEwLlIpDQoNCkZvciBSS0hTLCBoIGlzIGJhbmR3aWR0aCBwYXJhbWV0ZXIuIE5PVCBoZXJpdGliaWxpdHkuICANCg0KQmFzZWQgb24gW0NhbXBvcyBldCBhbC4sIDIwMTBdKGh0dHBzOi8vd3d3LmNhbWJyaWRnZS5vcmcvY29yZS9qb3VybmFscy9nZW5ldGljcy1yZXNlYXJjaC9hcnRpY2xlL3NlbWlwYXJhbWV0cmljLWdlbm9taWNlbmFibGVkLXByZWRpY3Rpb24tb2YtZ2VuZXRpYy12YWx1ZXMtdXNpbmctcmVwcm9kdWNpbmcta2VybmVsLWhpbGJlcnQtc3BhY2VzLW1ldGhvZHMvMkI4MjM5MTZDRjRENzZGQUU2QkM4MTQ1NUZENzNBQkIvY29yZS1yZWFkZXIpLCAgICANCg0KPiAiQXMgc2hvd24gaW4gRmlnLiA0LCB0aGUgdmFsdWUgb2YgdGhlIGJhbmR3aWR0aCBwYXJhbWV0ZXIgdGhhdCBnYXZlIHRoZSBiZXN0IHByZWRpY3RpdmUgYWJpbGl0eSB3YXMgaW4gdGhlIHJhbmdlICgyLDQpLCBleGNlcHQgZm9yIGVudmlyb25tZW50IEUyIGluIHdoaWNoIHZhbHVlcyBvZiBoIG5lYXIgb25lIHBlcmZvcm1lZCBzbGlnaHRseSBiZXR0ZXIuIiAgIA0KDQp3aGVyZSBmaWd1cmUgNCBpcyB0aGUgRXN0aW1hdGVkIE1TRSAobWVhbi1zcXVhcmVkIGVycm9yKSBiZXR3ZWVuIENWIHByZWRpY3Rpb25zICh5SGF0KSBhbmQgb2JzZXJ2YXRpb25zICh5KSB2ZXJzdXMgdmFsdWVzIG9mIHRoZSBiYW5kd2lkdGggcGFyYW1ldGVyLCBoLCBmcm9tIDAuMSAtIDEwLiBUaGUgYmFuZHdpZHRoIHBhcmFtZXRlciB3aXRoIHRoZSBsb3dlc3QgTVNFIGdhdmUgdGhlIGJlc3QgcHJlZGljdGl2ZSBhYmlsaXR5LiAgIA0KU28gYXMgd2UgbG9vayB0aGUgdGhlIGBQTVNFYCBwbG90IGFnYWluc3QgYGhgLCB3ZSBzaG91bGQgc2VsZWN0IHRoZSBiYW5kd2lkdGggcGFyYW1ldGVyIGBoYCB0aGF0IGdpdmVzIHRoZSBsb3dlc3QgTVNFIGZvciBlYWNoIG9mIHlvdXIgcGhlbm8gZGF0YXNldHMuICAgDQoNCkkgbmVlZCB0byBkZXRlcm1pbmUgd2hhdCBpcyB0aGUgb3B0aW1hbCBiYW5kd2lkdGggcGFyYW1ldGVyIHRvIHVzZSBmb3IgZWFjaCBkYXRhc2V0LiBbQ2FtcG9zIGV0IGFsLiwgMjAxMF0oaHR0cHM6Ly93d3cuY2FtYnJpZGdlLm9yZy9jb3JlL2pvdXJuYWxzL2dlbmV0aWNzLXJlc2VhcmNoL2FydGljbGUvc2VtaXBhcmFtZXRyaWMtZ2Vub21pY2VuYWJsZWQtcHJlZGljdGlvbi1vZi1nZW5ldGljLXZhbHVlcy11c2luZy1yZXByb2R1Y2luZy1rZXJuZWwtaGlsYmVydC1zcGFjZXMtbWV0aG9kcy8yQjgyMzkxNkNGNEQ3NkZBRTZCQzgxNDU1RkQ3M0FCQi9jb3JlLXJlYWRlcikgc2hvdWxkIHNoZWQgc29tZSBsaWdodCBvbiBob3cgdG8gc2VsZWN0IHRoZSBoLiAgIA0KDQpgYGB7UiBldmFsPUZBTFNFfSANCiMjIERldGVybWluaW5nIGJhbmR3aWR0aCBwYXJhbWV0ZXIgZm9yIFJLSFMgLS0tLS0tLS0tLS0tLS0tLQ0Kc2V0d2QoIkQ6L1BIUyBHZW5vbWljIFNlbGVjdGlvbiBQcm9qZWN0L0RhdGEgQW5hbHlzaXMvR1MvMjAxODEyMTgiKQ0KIyBzZXR3ZCgiL1VzZXJzL3NoYW50ZWwvRGVza3RvcC8yMDE4MTIxOCBHUyBNYWMiKQ0KZm9sZHM9NSAgICAgICAgICAgICAgDQpydW49NQ0Kbkl0ZXIgPSAxMjAwMA0KYnVybkluPTIwMDANCg0KI0JhbmR3aWR0aCBwYXJhbWV0ZXJzIHRvIGNvbXBhcmUNCmg9YyguMDEsLjEsLjQsLjgsMS41LDMsNSkNCg0KI1lvdXIgRGF0YSB+fn5+fn5+fn5+fn5+fn5+DQpteVkgPSBteVljOyBteUdEID0gbXlHRGMNCnkgPSBteVlbLDJdIA0KWCA9IG15R0Q7IG48LW5yb3coWCk7IHA8LW5jb2woWCkNCkQ8LShhcy5tYXRyaXgoZGlzdChYLG1ldGhvZD0nZXVjbGlkZWFuJykpXjIpL3ANCg0KI1N0YXRpc3RpYyBQYXJhbWV0ZXJzDQpQTVNFPC1udW1lcmljKCkgOyBWQVJFPC1udW1lcmljKCk7IFZBUlU8LW51bWVyaWMoKSA7DQpwRDwtbnVtZXJpYygpOyBESUM8LW51bWVyaWMoKQ0KZm1MaXN0PC1saXN0KCkNCg0KI1JLSFMgQ2FsY3VsYXRpb24gfn5+fn5+fn4NCmZvcihpIGluIDE6bGVuZ3RoKGgpKXsNCiAgSzwtZXhwKC1oW2ldKkQpDQogIEVUQTwtbGlzdChsaXN0KEs9Syxtb2RlbD0nUktIUycpKQ0KICBwcmVmaXg8LSBwYXN0ZShoW2ldLCAiXyIsc2VwPSIiKQ0KICBmbTwtQkdMUih5PXlOQSxFVEE9RVRBLCBuSXRlcj1uSXRlcixidXJuSW49YnVybkluLHNhdmVBdD1wcmVmaXgpDQogIGZtTGlzdFtbaV1dPC1mbQ0KICBQTVNFW2ldPC1tZWFuKCh5W3RzdF0tZm0keUhhdFt0c3RdKV4yKQ0KICBWQVJFW2ldPC1mbSR2YXJFDQogIFZBUlVbaV08LWZtJEVUQVtbMV1dJHZhclUNCiAgRElDW2ldPC1mbSRmaXQkRElDDQogIHBEW2ldPC1mbSRmaXQkcEQNCn0NCg0KI1Bsb3RzIGZvciBlYWNoIGggfn5+fn5+fg0KUjI8LTEtUE1TRS9tZWFuKCh5W3RzdF0tbWVhbih5Wy10c3RdKSleMikNCg0KcGxvdChWQVJFfmgseGxhYj0iQmFuZHdpZHRoIiwgDQogICAgIHlsYWI9IlJlc2lkdWFsIFZhcmlhbmNlIix0eXBlPSJvIixjb2w9NCkNCnBsb3QoSShWQVJFL1ZBUlUpfmgseGxhYj0iQmFuZHdpZHRoIiwNCiAgICAgeWxhYj0idmFyaWFuY2UgcmF0aW8gKG5vaXNlL3NpZ25hbCkiLHR5cGU9Im8iLGNvbD00KQ0KcGxvdChwRH5oLHhsYWI9IkJhbmR3aWR0aCIsIHlsYWI9InBEIix0eXBlPSJvIixjb2w9MikNCnBsb3QoRElDfmgseGxhYj0iQmFuZHdpZHRoIiwgeWxhYj0iRElDIix0eXBlPSJvIixjb2w9MikNCnBsb3QoUjJ+aCx4bGFiPSJCYW5kd2lkdGgiLCB5bGFiPSJSLXNxdWFyZWQiLHR5cGU9Im8iLGNvbD0yKQ0KcGxvdChQTVNFfmgseGxhYj0iQmFuZHdpZHRoIiwgeWxhYj0iUE1TRSIsdHlwZT0ibyIsY29sPTIpDQoNCmBgYA0KDQpBbHNvLCBJIGp1c3QgZm91bmQgYW5vdGhlciBhcnRpY2xlIGRpc2N1c3Npbmcgc2VsZWN0aW9uIG9mIHRoZSBiYW5kd2lkdGggcGFyYW1ldGVyIFtQ6XJlei1FbGl6YWxkZSBldCBhbC4sIDIwMTVdKGh0dHBzOi8vbGluay5zcHJpbmdlci5jb20vYXJ0aWNsZS8xMC4xMDA3L3MxMzI1My0wMTUtMDIyOS15KS4gTmVlZCB0byByZXZpZXcsIGJ1dCBmb3Igbm93LCBJIHRoaW5rIHdlIGhhdmUgYSBnb29kIHN0YXJ0IGZvciBoLg0KDQoNCg==

h: Bandwidth Parameter